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It has been reported recently that the equipartition theorem is violated in molecular dynamics 
simulations with periodic boundary condition [Shirts et al, J. Chem. Phys. 125 164102 (2006)]. 
This effect is associated with the conservation of the center of mass momentum. Here, we propose 
a fluctuating center of mass molecular dynamics approach (FCMMD) to solve this problem. Using 
the analogy to a system exchanging momentum with its surroundings, we work out -and validate 
via simulations- an expression for the rate at which fluctuations shall be added to the system. The 
restoration of equipartition within the FCMMD is then shown both at equilibrium as well as beyond 
equilibrium in the linear response regime. 



I. INTRODUCTION 

The equipartition theorem states that the total ki- 
netic energy of a classical system in canonical ensem- 
ble is equally distributed among all degrees of freedom 
and that the average kinetic energy associated with the 
translational motion of a particle is given by (jP /{2m)) = 
dk^T/2. Here, p and m are the momentum and mass of 
the particle, ks is the Boltzmann factor, T denotes the 
temperature and d is the spatial dimension. This relation 
serves to control the temperature in molecular dynamics 
(MD) simulations by adjusting the kinetic energy of the 
system [ij. 

It has been shown recently that in MD simulations with 
the periodic boundary condition (PBC) the equipartition 
theorem is violated [2|- This effect is attributed to the 
conservation of the center of mass (or, equivalently, to- 
tal) momentum, Pcm, due to the PBC. This additional 
constant of motion restricts the simulation trajectories 
to only a subset of the phase space and leads to a dif- 
ference between the time- and ensemble-averages [1, 0] ■ 
It can be shown that, in canonical ensemble MD simula- 
tions with PBC, the average kinetic energy of a particle 
of mass m obeys Q 
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2m 



dk^T 
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where 



stands for statistical average and A/t, 



J2i=i ""^i total mass of the system [N is the to- 

tal number of particles). Eq. ^ shows that the viola- 
tion of the equipartition theorem can be safely ignored 
if the mass of a single particle is negligible compared 
to the total mass of the system. Also, if all the parti- 
cles have the same mass, m/Mtot&i = l/N the effect is 
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negligible for most simulation cases. This is also in line 
with the general observation that differences between to 
molecular dynamics ensembles @ and other thermody- 
namic ensembles become less significant as the system 
size grows and eventually vanishes in the thermodynamic 
limit {N ^oo) 0. 

However, beyond iso-particle systems, the same prob- 
lem arises principally also in hetero-particle systems, 
where for instace a small number of massive particles are 
surrounded by a large number of light particles. This in- 
cludes explicit solvent MD simulations of transport prop- 
erties of colloids and nanoparticles in the dilute limit and 
generally all multi-component mixtures. In these cases, 
the violation of equipartitioning may not always be toler- 
able. For an estimate, let us consider a system consisting 
of A^i light particles of mass mo and massive parti- 
cles of mass amQ. This yields Aftotai = (-^i + Q;iVh)TOo 
and thus one obtains for the average kinetic energy of 
the heavy particle, {E^"') = [1 - a/{Ni + aN^)]dk^Tl2. 
Obviously, in the limit that the total mass of the heavy 
particles is large compared to the total mass of the 
light particles (a > iVi/iVh), this reduces to {E^™) w 
(1 — 1/N't,)dk-B,T /2. Thus, the average kinetic energy of 
a single heavy particle approaches zero with increasing 
mass. 

The present paper addresses this point. We propose 
a method to restore the equipartition theorem in molec- 
ular dynamics ensembles containating components with 
different masses. For this purpose, we first use computer 
simulations of a massive tracer particle in an ambient liq- 
uid and provide evidence for the idea that, as noted in 
@j 1 the main cause of the problem is the conservation 
of the center of mass momentum. Based on this un- 
derstanding, wc propose a molecular dynamics method 
which allows for the fluctuations of the center of mass 
momentum, Pcm (note that, in equilibrium, (5JPcm = -Pcm, 
since (Pcm) = 0). For a reliable implementation of the 
method, we also determine the rate at which fluctua- 
tions shall be added to the system, Eq. ([5]). The va- 
lidity of this expression is conflrmed via computer sim- 
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ulations of systems with rigid walls, which do not have 
the artificial constant motion, i-'cm = 0. Finally, we test 
the method showing that it does restore equipartitioning 
both in equilibrium and beyond equilibrium in the linear 
response regime. 



II. VIOLATION OF EQUIPARTITION FOR A 
MASSIVE TRACER 



As mentioned above, a strong violation of the equipar- 
tition theorem is expected for the case of a massive tracer 
particle in a liquid environment. In order to demon- 
strate this property, we perform MD simulations of a 
generic 80:20 binary mixture of Lennard- Jones parti- 
cles (types A and B) 0, A and B particles in- 
teract via ULjir) = 4e„^[(da^/r)i2 - (dap/rf], with 
a, /? = A,B, EAB l-5eAA, ebb = O.Seaa, ^ab = O.S^aa, 
^BB = O.SSdAA; and toa — ttib- The potential is trun- 
cated at twice the minimum position of the LJ potential, 
^c,Q/3 = 2.2A5da/3- The parameters caa, ^aa and toa 
define the units of energy, length an d mass. The unit of 
time is given by tlj = cJaa -vA^a/^aa ■ The system density 
is kept constant at the value of 1.2 and temperature at 
T = 1 for all simulations whose results are reported here. 
Depending on the case studied, linear dimension and to- 
tal particle number arc in the range of L e [4.7, 203] and 
N G [125,10684]. Equations of motion are integrated 
using the velocity- Verlct algorithm with a discrete time 
step of (it = 0.005. 

The results presented here are expected to be largely 
model independent and hence general. The choice of the 
above model is purely historical and is motivated by the 
fact that we have been using it to study a number of 
problems in the context of the physics of glasses p^| - [l3j . 
We indeed encountered the present problem of the viola- 
tion of the equipartition theorem as we inserted massive 
tracer particles into our model to study the concept of 
effective temperature [13, [l^ . 

With the exception of one, the mass of all particles is 
set to unity. One of the particles (of type B) is taken to 
be the massive tracer. The mass of this particle is then 
varied and its kinetic energy is monitored. Simulation 
results are averaged over 40 independent runs. In or- 
der to investigate the possible role of the thermostat, all 
the simulations are performed both for the Nose-Hoover 
(N-H) [H [i3 and the Andersen ^ thermostats. In 
equilibrium simulations, all components of particle veloc- 
ities are coupled to the thermostat. We also extend the 
present analysis to a non-equilibrium steady state situ- 
ation by imposing a linear shear flow via the SLLOD- 
algorithm combined with the Lees-Edwards boundary 
condition (LEBC) [1^]. In this case, coupling to the 
thermostat is done only for the velocity component in 
the direction perpendicular to the shear plane (vorticity 
direction). By doing so, we avoid problems related to the 
flow-induced bias on the kinetic energy, when regulating 
the system temperature [20| . 
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FIG. 1. (a) The average kinetic energy of a massive parti- 
cle, normalized by dk^T /2, is shown versus its mass for two 
popular thermostating methods, both in equilibrium and un- 
der steady shear (with uniform shear rate of 7 = 10~^). The 
solid line gives the analytic prediction given by Eq. ([T]). Re- 
sults for a system containing two planar walls are also shown 
for equilibrium simulations. While equipartition is violated 
in the simulations with PBC, it is well satisfied in the simula- 
tions with walls, (b) Total momentum versus time for all the 
simulations reported in the panel (a). The inset shows Pcm- 
distribution for the case of simulation with walls. It perfectly 
obeys the expected Gaussian distribution with zero mean and 
variance MtotaifceT (The black dots are simulation results for 
Mtotai = 125 and the continues line is the corresponding the- 
oretical expectation). 



Results obtained via these simulations are depicted in 
Fig. ([T]). As seen in this plot, the violation of the equipar- 
tition theorem occurs in perfect agreement with the the- 
oretical predictions of Eq. ([1]), independent of the specific 
thermostat. 

Next we provide evidence from simulation that the 
deep reason for the violation of the equipartition theo- 
rem is indeed the conservation of the total momentum 
0, d] . For this purpose, we have designed a simulation 
setup where Pcm is not conserved. This is achieved by 
introducing two planar walls separated by a distance 
along the z-direction, while PBC is used along the x and 
y directions. The walls arc made of particles with the 
same size and structure as the liquid particles so that 
liquid-wall interactions induce fluctuations of Pcm along 
all spatial directions. Results of these simulations are 
also shown in Fig. ([1]), demonstrating that, as expected, 
the equipartition theorem is valid in systems with walls. 



III. FLUCTUATING CENTER OF MASS MD 
(FCMMD) 

The above results suggest that a possibility to re- 
store the equipartition theorem is to introduce walls with 
roughness on the particle scale. However, in studies fo- 
cusing on bulk properties, walls are undcsircd since they 
in general influence the system properties unless very 
large wall-to-wall separations are used (see, e.g., j2l| - [23| 
and references therein). Thus, it is desirable to intro- 
duce a method which uses PBC, while at the same time 
allowing for fluctuations of the total momentum. Such a 
method is proposed here. Our approach is quite simple 
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and is motivated by the fact that, in a system exchang- 
ing momentum with its environment, the center of mass 
momentum is a fluctuating quantity. 

Motivated by this idea, we perform the following two 
steps: (i) Draw a value for Pcm and (ii) distribute it 
among particles. These steps are carried over, repeat- 
edly, during the simulation. In order to have the canon- 
ical sampling of the phase space, the total momentum 
assigned in step (i) should assume a distribution proba- 
bility coinciding with the canonical distribution function 
for P™: 



N 



27rfcBTMtotal 



exp 



2fcBTA/total 



(2) 



where /(x) is the probability of the micro-state Xi and 
S is the Dirac delta-function. Similar to the discussion 
in Ref. [l[ for choosing kinetic energy from the canoni- 
cal distribution, one has a certain flexibility in choosing 
the sampling rate. Here, we provide a physical criterion 
to estimate the time scale of the Pcm-fiuctuations. We 
build our analysis upon the fact that fluctuations of Pan 
are caused by the exchange of momentum with the sur- 
rounding medium. 

For a system with Pcm{t = 0) = 0, it follows from the 
above considerations that, due to interactions with the 
surrounding medium, Pcm will not remain zero but in- 
creases with time. On the other hand, too large a value 
of Pcm will decay due to the same interactions. Collisions 
with the surrounding medium thus provide a source of 
stochastic noise and, at the same time, give rise to vis- 
cous friction. This is very similar to the fluctuations of 
the velocity of a Brownian particle in a fluid. The proba- 
bility distribution of these fluctuations is obtained as the 
solution of a Fokker-Planck equation subjected to the 
potential Q , 



af(Pcm,0 

dt 



MV|p„-(fV|p,»+i5(V|p„)'f, (3) 



where 0(Pcn,) = -P'i/ {2M^,^,^) and V|p„ ^ 
{d/dPcm,x,d/dPcm,y,d/dPcmz,z)- Thc mobility, fi, and 
diffusion constant, D, obey the Einstein relation, fi = 
D/{k^T). Given P^'^ at time t' , the conditional probabil- 
ity distribution at a time t > t' is [131 > 



f(p,„,t|p,;,t') 



1 



(Pc„-P,;cxp[-7(t-t')])' 



X exp 



(4) 



where 7 = /i/Mtotai and cr^(t,t') = dhsT Mtotaii^ - 
exp(— 27(t — t'))]. It is seen from Eq. that p reaches 
the expected Maxwell distribution, Eq. ([2]), in the limit 



of long times, t — t' ^ I/7. The characteristic time for 
reaching the equilibrium distribution of center of mass 
fluctuations is thus obtained from r = I/7, 



D 



(5) 



This expression is not fully satisfactory as it contains 
an important unknown parameter, D. We therefore at- 
tempt at an estimate of r from a microscopic consid- 
eration. For this purpose, we again recall that, start- 
ing with Pcm(i = 0) = collisions with the surrounding 
medium will lead to (Pcm) = (iMtotai^B?" within a time 
of the order of r. For simplicity, we assume here that 
Pcm is the sum of Ns statistically independent elemen- 
tary momentum fluctuations, 6pi, resulting from the col- 
lisions between fluid particles with the system's bound- 
ary, P,^ = ^f' Sp,. This yields (P^) = N, {6p^) and 
thus Ns {Sp^'^ = dAItotaik^T . Thc time scale r is en- 
coded in the number of elementary collisions Ng. To 
see this, we first note that momentum exchange occurs 
within a "skin" - which runs parallel to the boundary 
- of thickness equal to mean free path, ^froo- On aver- 
age, 1/6*^ of these particles in the skin layer move along 
the perpendicular direction toward the boundary and will 
undergo a collision within a time of 5t ~ htce / {\v = 
IhccI \fk^Fjra where (|f±|) is the average thermal veloc- 
ity in the direction normal to the boundary. The total 
number of collisions within a time of r is thus obtained 
as ^ p/frec^/S X Tjbt ^ pATy/k^rJra/&. To ar- 
rive at a closed expression for r, the magnitude of the 
typical momentum exchange per collision is estimated: 
5p ^ 2to(|i)j^|) = 2\/mk-e,T. Combining the above two 
expressions for and using this last relation for 5p, one 
finally finds 



C 



T = 



p{mkBTy/^ A 



(6) 



where C is a constant prefactor. Equation ^ gives an es- 
timate for the characteristic time of the Pcm-fluctuations 
in a system exchanging momentum with its surroundings 
through a boundary (interface) of surface area A. 

In order to test this result, we have performed a se- 
ries of three dimensional MD simulations of the present 
binary LJ model confined between two parallel walls for 
different system sizes while keeping all other simulation 
parameters constant (e.g., T = 1, p = 1.2, m = 1). The 
characteristic time is measured by the auto-correlation 
time of fluctuations of Pcm. Note that, in these simula- 
tions, PBC is used along the x and y directions, so that 
no momentum fluctuations will originate from the corre- 
sponding boundaries. In other words, the relevant sur- 
face area. A, appearing in Eq. ^ corresponds to the sur- 
face area of the walls. To better highlight the dependence 
of T on Mtotai and A, we studied two different geometries 
leading to qualitatively different results for r in terms of 
the total mass. In the flrst series of simulations, the sys- 
tem was a cube with length L so that Mtotai = pL^ and 
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FIG. 2. The time scale, r, for Pcm-fluctuations versus to- 
tal mass in a system confined by two planar walls, r is 
determined from the decay of the autocorrelation function, 
(fern (t) fern (0)). Results are shown for two different geome- 
tries. In case I, the simulation box is a cube and the variation 
of Aftotai is accompanied by a corresponding change of the 
surface area of the walls. In case II, the surface area of the 
walls is kept constant but only the distance of the walls is 
varied. Using Eq. ((6|, one thus expects r oc M^l^^-^ in case I 
but T oc Mtotai in case II (solid lines). 
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FIG. 3. Restoring the equipartition theorem via the proposed 
fluctuating center of mass molecular dynamics (FCMMD) 
method. The average kinetic energy of a massive particle, 
normalized by d k^T 12, is shown versus its mass both for equi- 
librium simulations as well as under steady shear in the linear 
response regime (uniform shear rate of 7 = 10"''). For com- 
parison, we also plot results of MD simulations without Pcm- 
fluctuations and the corresponding theoretical curve, Eq. ([1]). 



IV. CONCLUSION 



A = L'^ = (Mtotai/p)^/^ and thus t a M^/^i (case I). In 
the second series of simulations, we only varied the wall- 
to-wall separation, i^, while keeping the surface area of 
the walls constant. This gives r ex Mtotai (case II). As 
shown in Fig. ([2]), results on the characteristic time of 
momentum exchange obtained for these two sets of sim- 
ulations clearly confirm the validity of Eq. ([B]) with a 
constant of proportionality of C ~ C(l). 

In our scheme, we update P^a using a random walk 
sampling f(Pcm) with time scale of order r, Eq. ([6|). The 
question arises now as how to distribute a given Pcm 
among particles. In canonical ensemble, if the total mo- 
mentum of the system is J'cm, the conditional average 
for the momentum of a particle is equal to (p)|f>„ = 
mPcm/Mtotai- Wc, therefore, propose that the imposed 
momentum change be divided among particles propor- 
tional to their individual masses. In analogy to the sys- 
tem with walls, the momentum change is applied to each 
particle once at a time. The order of particles is chosen 
randomly. After adding momentum to each particle, rel- 
ative velocities of all particles with respect to the center 
of mass are rescaled. This last operation does not modify 
the center of mass momentum but allows to restore the 
kinetic energy exactly to the value before updating Pcm- 

Results obtained from these simulations are shown in 
Fig. ([3]). As shown in this figure, the proposed approach 
is able to restore the equipartition theorem both in equi- 
librium simulations as well as in a system beyond equi- 
librium in the linear response regime. 



In this work, we propose a modification of the molec- 
ular dynamics method with periodic boundary condition 
to restore the equipartition theorem. The method is 
based on introducing fiuctuations of the center of mass 
momentum. The issue of a proper rate at which Pcm 
fluctuations are imposed to the system is also addressed 
and validated against simulations. It is shown that the 
method restores equipartition both at equilibrium and 
under steady shear in the linear response regime. This 
latter flnding is of crucial importance for studies, which 
focus on a violation of the equipartition due to non-linear 
off-equilibrium effects [13, Hg • It is noteworthy that the 
violation of equipartition does not exclusively occur in 
MD simulations. As an example, it has also been ob- 
served in the fluctuating lattice Boltzmann method where 
equipartitioning is important at all length scales [25j . 
The relevance of the present work is thus not restricted 
to MD simulations but may also provide guidance for 
restoring equipartitioning and hence, a correct thermo- 
stat method, in other mesoscale simulations [26l. [27|. 

V. ACKNOWLEDGMENTS 

Nima H. Siboni gratefully acknowledges the finan- 
cial support from the Deutsche Forschungsgemeinschaft 
(German Research Foundation) through grant GSC 111. 
ICAMS gratefully acknowledges funding from its indus- 
trial sponsors, the state of North-Rhine Westphalia and 
the European Commission in the framework of the Eu- 
ropean Regional Development Fund (ERDF). 



5 



G. Bussi. D. Donadio. and M. Parrinell o. fThe Journal 
of Chemical Physics, 126, 014101 (2007)[ 



[16] S. Nose, The Journal of Chemical Physics. 81. 511 (1984) 



R. B. Shirts. S. R. Burt, and A. M. Johnso n. The Journal 
of Chemical Physics, 125 164102 (2006)[ 



[17] W. G. Hoover, Phys. Rev. A, 31, 1695 (T98F 



J. R. Ray and H. Zhang, [Phys. Rev. E, 59, 4781 (1999) 
M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Mar- 



M. J. Uline, D. W. Siderius, and D. S. Corti, The Journal 



[18] H. C. Andersen, J. Chem. Phys., 72, 2384 (1980). 
[19] D. J. Evans and G. P. Morriss, Statistical Mechanics 
of Non Equilibrium Liquids (Academic Press, London, 
1990) 



tvna. The Journal of Chemical Phvsics. 115. 1678 f2001' . [20] D. J. Evans and G. P. Morriss, Phys. Rev. Lett., 56, 2172 



(1986) 



of Chemical Physics, 128, 124301 (2008) 
J. J. Erpenbeck and W. W. Wood, Statistical Mechanics 
Part B: Time dependent processes, edited by B. J. Berne 
(Plenum Press, 1977) ISBN 0306335050 



J. L. Lebowitz. J . K. Percus, and L. Verlet, Phys. Rev 



153, 250 (1967) 



W. Kob and H. C. Andersen, Phys. Rev. Lett., 73, 1376 
(1994). 

W. Kob and H. C. Andersen, Phys. Rev. E, 51, 4626 
fl995l. 



[21] F. Varnik, J. Baschnagel, and K. Binder, Phy. Rev. E, 

65, 021507 (2002). 
[22] F. Varnik, J. Baschnagel, and K. Binder, Eur. Phys. J. 

E, 8, 175 (2002). 
[23] J. Baschnagel and F. Varnik, J. Phys.: Condens. Matter, 

17, R851 (2005). 
[24] H. Risken, The Fokker-Planck Equation, 2nd ed. 

(Springer, 1989). 
[25] S. T. T. Omia, C. Denniston, M. Karttunen. and T. Ala- 



Nissila. 



F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys. 



(2011) 



120, 2788 (2004). 

F. Varnik, J. Chem. Phys., 125, 164514 (2006). 

F. Varnik and O. Henrich. Phy. Rev. B. 73, 174209 



The Journal of Chemical Physics, 134, 064902 



[26] M. Gross, R. Adhikari, M. E. G ates, and F. Varnik, 



Phys. Rev. E, 82, 056714 (2010) 



[27] M. Gross, M. Gates, F. Va rnik, and R. Adhikari, J. Stat 



(2006). 



F. Varnik and D. Raabe, phys. Rev. E, 77, 011504 (2008)| . 

P. SoUich, Phys. Rev. E, 58, 738 (1998). 

L. Berthier and J.-L. Barrat, J. Chem. Phys., 116, 6228 

(2002). 



Mech., 03, P03030 (2011) 



